%% bdary_face (for Neumann and Robin)
function bdf = bdary_face(mesh, index)
if isfield(mesh, "bdary")
    bdf = mesh.bdary(mesh.bdary(:,5) == index,:);
else
    bdf = mesh.bdf(mesh.bdf(:,5) == index,:);
end
end